subroutine CalcFacesVariables()
use Data_cfoam
use Data_Mesh
use DataTime, only: dt, nstep, ttime, dt2
use DataPerfectGas, only: cv, Gamma, TemperaturePRho
use DataBoundarySource, only: BoundSourcePresent, iBoundSurf, CalcBoundSource
use InletBoundaryConditions, only: InletBCon
implicit none
include 'mpif.h'
integer :: i,j,k, l,m, j1, iFace
integer :: ndeltap,bmeshs,bmeshf
integer :: iProc,statusp(mpi_status_size) 
real(8) :: dtminx
real *8 uf(3), ucf(3)  ! local velocities' basis
real *8 ucsf(3) ! old value of ucf ! s = staroe!
real(8) :: alfa, NoSlipCoef
real *8     cZ(3) ! coord of cell centers
real(8) :: rr

real(8) :: R(3), rMod,c, time, Q,dQdT,p0, a
real(8) :: U0(3), Pres0, Temp0, Rho0
real(8) :: finlet1

real(8) :: r1(3), pt1, pr(3), gradP(3), ut1(3),ur(3,3), divU, Err(3), pt,ut(3), dv, ErrPt
real(8) :: x,y,z, Umean, Lr, Px, CharN(3), s1, s2

! ALCF: requires verification
!$omp parallel do private( i,j,l,m,CharN,s1,s2 )
do k=1,nInternalFaces
    i = meshFaceOwner(k)
    j = meshFaceNeighbour(k)
    l = opposingOwnerFaceLabel(k)
    m = opposingNeighbourFaceLabel(k)
            
    if(NormalChars)then
        CharN(1:3) = meshSf(1:3,k)
    else
        CharN(1:3) = meshCellCenters(1:3,i) - meshCellCenters(1:3,j)
    end if
    
    call usdefn(Usnew(1:3,k),Us(1:3,k),Us(1:3,l),Us(1:3,m),U(1:3,i),Unew(1:3,i),U(1:3,j),Unew(1:3,j) &
           ,Psnew(k),    Ps(k),    Ps(l),    Ps(m),    p(i),    pnew(i),    p(j),    pnew(j)         &
           ,Rhosnew(k),  Rhos(k),  Rhos(l),  Rhos(m),  rho(i),  rhonew(i),  rho(j),  rhonew(j)       &
           ,Db(k),Df(k),CharN,dt,dtminx,epsilonSmoothed)
    
    ! ALCF modifed, Phil, 2012-02-28
    s1 = ( gam - 1d0 ) * Rhosnew( k )
    s2 = Usnew( 1,k ) * Usnew( 1,k ) + Usnew( 2,k ) * Usnew( 2,k ) + Usnew( 3,k ) * Usnew( 3,k )
    Esnew(k) = Psnew(k) / s1 + 0.5d0 * s2
end do



    include "netInterface.h"
    do j=1,nBoundaryFields
        if (boundaryMeshName(j)(1:4).eq.'proc') cycle ! goto 10
        bmeshs = 1 + boundaryMeshStart(j) ! Phil, 2012-02-28
        bmeshf = boundaryMeshStart(j) + boundaryMeshSize(j)
        
        if (btype(j).eq.6) then ! wall, full-Slip
            do k = bmeshs,bmeshf
                i=meshFaceOwner(k)
                
                Usnew(1:3,k) = U(1:3,i)
                Psnew(k) = P(i)
                Rhosnew(k) = Rho(i)
                Esnew(k) = Enew(i)
            end do
        end if
        
        if (btype(j).eq.4) then ! wall, full-Slip
            do k = bmeshs,bmeshf ! 1+boundaryMeshStart(j),boundaryMeshStart(j)+boundaryMeshSize(j)  ! Phil, 2012-02-28
                ! SK changed this to a standard wall bcs..
                i=meshFaceOwner(k)
                l=opposingOwnerFaceLabel(k)
                ! SK changed here
!                    utmp(6,100)=epsilon
!                    ntimemax=20000
!                    if(nstep.lt.ntimemax) then
!                        propor=(nstep-1.0)/(ntimemax-1.0)
!                        utmp(6,100)=1.0d0-propor+propor*epsilon
!                    end if
                NoSlipCoef = 1.0d0
                call usdef3(Usnew(1:3,k),Us(1:3,k),Us(1:3,l),U(1:3,i),Unew(1:3,i) &
                       ,Psnew(k),    Ps(k),    Ps(l),     p(i),    pnew(i)    &
                       ,Rhosnew(k),  Rhos(k),  Rhos(l),   rho(i),  rhonew(i)  &
                !	           ,tsnew(k),    ts(k),    ts(l),     t(i),    tnew(i)  &
                       ,Db(k),meshSf(1:3,k),dt,dtminx,epsilonSmoothed,NoSlipCoef)
                Esnew(k)=Psnew(k)*((gam-1.0d0)*Rhosnew(k))**(-1)+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)*0.5d0 ! Phil, 2012-02-28
                ! Esnew(k)=Psnew(k)/((gam-1.0d0)*Rhosnew(k))+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)/2.0d0 
            end do 
        end if
        if (btype(j).eq.1) then ! Outlet
            do k = bmeshs,bmeshf ! 1+boundaryMeshStart(j),boundaryMeshStart(j)+boundaryMeshSize(j)  ! Phil, 2012-02-28
                i=meshFaceOwner(k)
                l=opposingOwnerFaceLabel(k)
                ! SK changed here
!                utmp(4,100)=0.0
!!____________________________________________________________________________________________!
                if(ViscPartOnly)then
                    x = meshFaceCenters(1,k)
                    y = meshFaceCenters(2,k) !- 0.5d0
                    z = meshFaceCenters(3,k)
                    
                    Umean = 0.01
                    Lr = 0.005
                    
                    Px = - 12d0 * 1.5e-5 / Lr**2 * Umean
                    
                    utmp(5,j) = Pfon + Px * x
                    utmp(1,j) = -0.5d0 / 1.5e-5 * Px * z*(Lr-z)
                    utmp(4,j) = Tfon
                end if
!!                r(1:3) = meshFaceCenters(1:3,k)
!!                rMod = sqrt(sum(r(1:3)*r(1:3)))
!!                c = sqrt(1.4*PFon/RoFon)
!!                time = ttime + dt2
!!                call CalcQ3D(r,time,Q,dQdt,p0,U0)                
!!
!!                Rho0 = RoFon + (P0-PFon)/c**2
!!                Temp0 = p0 / (Rho0 * cv * (Gamma-1d0))
!!                Pres0 = p0
!!____________________________________________________________________________________________!
                if(TotalPresTempGiven)then
                    call CalcOutletBoundCond_TotalPresTempGiven(Usnew(1:3,k),Psnew(k),Rhosnew(k), &
                                    pnew(i),Rhonew(i),utmp(1,j))
!                    Usnew(1,k) = utmp(1,j)
!                    Usnew(2:3,k) = 0d0
!                    Psnew(k) = pnew(i)
!                    Rhosnew(k) = Rhonew(i)
                else
                    U0(1) = utmp(1,j) ! UFON
                    U0(2) = utmp(2,j) ! VFON
                    U0(3) = utmp(3,j) ! WFON
                    Temp0 = utmp(4,j)
                    Pres0 = utmp(5,j)
!***********************************************************************************!
!                    r(1:3) = meshFaceCenters(1:3,k)
!!                    call CalcQ2D_Hankel(r,ttime+dt2,Pres0,U0,Rho0)
!                    call CalcQ3D_ConvectedMonopole(r,ttime+dt2,Pres0,U0,Rho0)
!
!                    Temp0 = Pres0 / (Rho0 * cv * (Gamma-1d0))
!***********************************************************************************!
                   ! if(sum(U0*meshSf(1:3,k)).lt.0)then  ! always outlet
                    !    finlet1 = 1d0
                    !else
                        finlet1 = 0d0
                    !end if
                    call usdef1(Usnew(1:3,k),Us(1:3,k),Us(1:3,l),U(1:3,i),Unew(1:3,i) &
                       ,Psnew(k),    Ps(k),    Ps(l),     p(i),    pnew(i)    &
                !	           ,tsnew(k),    ts(k),    ts(l),     t(i),    tnew(i)  &
                       ,Rhosnew(k),  Rhos(k),  Rhos(l),   rho(i),  rhonew(i)  &
                       ,Db(k),meshSf(1:3,k),dt,dtminx,epsilonSmoothed,Pres0,Temp0,U0(1:3),finlet1) !,0.d0) !utmp(4,100)) !,alfa) 
                end if
                Esnew(k)=Psnew(k)*((gam-1.0d0)*Rhosnew(k))**(-1)+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)*0.5d0 ! Phil, 2012-02-28
                ! Esnew(k)=Psnew(k)/((gam-1.0d0)*Rhosnew(k))+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)/2.0d0  
            end do
        end if
        if (btype(j).eq.2) then ! Inlet
        !	call interp(p1,np1,ttime,uin)
            do k = bmeshs,bmeshf ! 1+boundaryMeshStart(j),boundaryMeshStart(j)+boundaryMeshSize(j)  ! Phil, 2012-02-28
                i=meshFaceOwner(k)
                l=opposingOwnerFaceLabel(k)
                cZ(1:3) = meshFaceCenters(1:3,k)
                if(ViscPartOnly)then
                    x = meshFaceCenters(1,k)
                    y = meshFaceCenters(2,k) !- 0.5d0
                    z = meshFaceCenters(3,k)
                    
                    Umean = 0.01
                    Lr = 0.005
                    Px = - 12d0 * 1.5e-5 / Lr**2 * Umean
                    
                    utmp(5,j) = Pfon + Px * x
                    utmp(1,j) = -0.5d0 / 1.5e-5 * Px * z*(Lr-z)
                    utmp(4,j) = Tfon
!                    Pin = 1.0001d5
!                    Pout = 1d5
!
!                    x = meshFaceCenters(1,k)
!                    y = meshFaceCenters(2,k) !- 0.5d0
!                    utmp(5,j) = Pin + (Pout-Pin) * x
!!                    a = 0.5 * (Pout-Pin) !/ nu
!                    utmp(1,j) = 0.5 * (Pout-Pin) *(y**2 - 0.25d0)
!                    !Us(1,i) = 0.5 * (Pout-Pin) * (y-1d0) * y / nu
!                    utmp(2,j) = 0d0
!                    utmp(3,j) = 0d0
!                    utmp(4,j) = Tfon
                end if
                U0(1) = utmp(1,j)
                U0(2) = utmp(2,j)
                U0(3) = utmp(3,j)
                Temp0 = utmp(4,j)
                Pres0 = utmp(5,j)
! SK we don't need this interpolation ?
!                if(PS90)then
                if(InletBCon)then
                    cZ(1:3) = meshFaceCenters(1:3,k)
                    rr=dsqrt(cZ(2)**2+cZ(3)**2)
!!         if (rr.le.0.7d0) &
!                    if (rr.le.0.025d0) &
!                        call intervt(xu,16,rr,utmp(1,100),utmp(3,100),utmp(2,100))
                        
!                    if (rr.le.0.025d0) then
                        call CalcNozzleBoundCond(cZ,U0,Temp0,Pres0)
!                    end if
                end if
!                utmp(4,100)=1.0d0
!                    utmp(6,100)=epsilon
                if(nstep.lt.nTimeStepAcceleration) then ! ntimemax
                    cZ(1:3) = meshFaceCenters(1:3,k)
                    U0(1) = U0(1)*propor+(1.-propor)*UFON
                    U0(2) = U0(2)*propor+(1.-propor)*VFON
                    U0(3) = U0(3)*propor+(1.-propor)*WFON
                    Pres0 = Pres0*propor+(1.-propor)*PFON
                    Temp0 = Temp0*propor+(1.-propor)*TFON
!                    utmp(1,100)=utmp(1,j)*propor+(1.-propor)*UFON
!                    utmp(2,100)=utmp(5,j)*propor+(1.-propor)*PFON
!                    utmp(3,100)=utmp(4,j)*propor+(1.-propor)*TFON
!!                    else
!!                        utmp(1,100)=utmp(1,100)*cos(0.001*(nstep-ntimemax))
                end if
!____________________________________________________________________________________________!
!                r(1:3) = meshFaceCenters(1:3,k)
!!                rMod = sqrt(sum(r(1:3)*r(1:3)))
!!                c = sqrt(1.4*PFon/RoFon)
!!                time = ttime + dt2
!                call CalcQ2D_Hankel(r,ttime+dt2,Pres0,U0,Rho0)
!!                call CalcQ3D(r,time,Q,dQdt,p0,U0)                
!!
!                Temp0 = Pres0 / (Rho0 * cv * (Gamma-1d0))
!                Pres0 = p0
!                Rho0 = RoFon + (P0-PFon)/c**2
!                Temp0 = p0 / (Rho0 * cv * (Gamma-1d0))
!!!                utmp(1,100) = U0(1)
!____________________________________________________________________________________________!
                if(TotalPresTempGiven)then
                    call CalcInletBoundCond_TotalPresTempGiven(Usnew(1:3,k),Unew(1:3,i), &
                                            Psnew(k),Rhosnew(k),meshSf(1:3,k),utmp(2,100),utmp(3,100),utmp(1,100))
                else
                   ! if(sum(U0*meshSf(1:3,k)).lt.0)then   ! always inlet
                        finlet1 = 1d0
                    ! else
                    !    finlet1 = 0d0
                    ! end if
                    call usdef1(Usnew(1:3,k),Us(1:3,k),Us(1:3,l),U(1:3,i),Unew(1:3,i) &
                            ,Psnew(k),    Ps(k),    Ps(l),     p(i),    pnew(i)    &
                            !	           ,tsnew(k),    ts(k),    ts(l),     t(i),    tnew(i)  &
                            ,Rhosnew(k),  Rhos(k),  Rhos(l),   rho(i),  rhonew(i)  &
                            ,Db(k),meshSf(1:3,k),dt,dtminx,epsilonSmoothed,Pres0,Temp0, & ! utmp(1,100),utmp(4,100)) !,alfa)  
                             U0(1:3),finlet1) ! 1d0) ! utmp(4,100)) !,alfa)  
                end if
!                call usdef1_TP(Usnew(1:3,k),Us(1:3,k),Us(1:3,l),U(1:3,i),Unew(1:3,i) &
!                            ,Psnew(k),    Ps(k),    Ps(l),     p(i),    pnew(i)    &
!                            !	           ,tsnew(k),    ts(k),    ts(l),     t(i),    tnew(i)  &
!                            ,Rhosnew(k),  Rhos(k),  Rhos(l),   rho(i),  rhonew(i)  &
!                            ,Db(k),meshSf(1:3,k),dt,dtminx,epsilonSmoothed,utmp(2,100),utmp(3,100),utmp(1,100),utmp(4,100)) !,alfa)  
                Esnew(k)=Psnew(k)*((gam-1.0d0)*Rhosnew(k))**(-1)+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)*0.5d0 ! Phil, 2012-02-28
                ! Esnew(k)=Psnew(k)/((gam-1.0d0)*Rhosnew(k))+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)/2.0d0 
            end do 
        end if
        if (btype(j).eq.3) then ! Symetry plane ! wall
            do k = bmeshs,bmeshf ! 1+boundaryMeshStart(j),boundaryMeshStart(j)+boundaryMeshSize(j)  ! Phil, 2012-02-28
                i=meshFaceOwner(k)
                l=opposingOwnerFaceLabel(k)
            ! SK changed here
                NoSlipCoef = 0.0d0
                call usdef3(Usnew(1:3,k),Us(1:3,k),Us(1:3,l),U(1:3,i),Unew(1:3,i) &
                           ,Psnew(k),    Ps(k),    Ps(l),     p(i),    pnew(i)    &
                           ,Rhosnew(k),  Rhos(k),  Rhos(l),   rho(i),  rhonew(i)  &
            !	           ,tsnew(k),    ts(k),    ts(l),     t(i),    tnew(i)  &
                           ,Db(k),meshSf(1:3,k),dt,dtminx,epsilonSmoothed,NoSlipCoef)
                Esnew(k)=Psnew(k)*((gam-1.0d0)*Rhosnew(k))**(-1)+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)*0.5d0 ! Phil, 2012-02-28
                ! Esnew(k)=Psnew(k)/((gam-1.0d0)*Rhosnew(k))+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)/2.0d0 
            end do 
        end if
        if (btype(j).eq.5) then ! Periodic
            ndeltap=boundaryMeshSize(j)/2
            do k=1+boundaryMeshStart(j),boundaryMeshStart(j)+ndeltap
                i=meshFaceOwner(k)
                j1=meshFaceOwner(k+ndeltap)
                l=opposingOwnerFaceLabel(k)
                m=opposingOwnerFaceLabel(k+ndeltap)
                Df(k)=Db(k+ndeltap)
                ! SK changed here
                call usdefn(Usnew(1:3,k),Us(1:3,k),Us(1:3,l),Us(1:3,m),U(1:3,i),Unew(1:3,i),U(1:3,j1),Unew(1:3,j1) &
                       ,Psnew(k),    Ps(k),    Ps(l),    Ps(m),    p(i),    pnew(i),    p(j1),    pnew(j1)     &
                       ,Rhosnew(k),  Rhos(k),  Rhos(l),  Rhos(m),  rho(i),  rhonew(i),  rho(j1),  rhonew(j1)     &
                !	           ,Tsnew(k),    Ts(k),    Ts(l),    Ts(m),    T(i),    Tnew(i),    T(j1),    Tnew(j1)     &
                       ,Db(k),Df(k),meshSf(1:3,k),dt,dtminx,epsilonSmoothed)
                Esnew(k)=Psnew(k)*((gam-1.0d0)*Rhosnew(k))**(-1)+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)*0.5d0 ! Phil, 2012-02-28
!                Esnew(k)=Psnew(k)/((gam-1.0d0)*Rhosnew(k))+(Usnew(1,k)**2+Usnew(2,k)**2+Usnew(3,k)**2)/2.0d0 
                Usnew(1:3,k+ndeltap)=Usnew(1:3,k)
                Psnew(k+ndeltap)=Psnew(k)
!                Tsnew(k+ndeltap)=Tsnew(k)
                Rhosnew(k+ndeltap)=Rhosnew(k)
                Esnew(k+ndeltap)=Esnew(k)
            end do
        end if
! 10         continue
    end do
    
    if(BoundSourcePresent)then
        j = iBoundSurf
        bmeshs = 1 + boundaryMeshStart(j) ! Phil, 2012-02-28
        bmeshf = boundaryMeshStart(j) + boundaryMeshSize(j)
        do iFace = bmeshs,bmeshf
            Usnew(3,iFace) = Usnew(3,iFace) + CalcBoundSource(meshFaceCenters(1:3,iFace),ttime+dt)
            a = CalcBoundSource(meshFaceCenters(1:3,iFace),ttime+dt)
            if(abs(a).gt.1d-5)then
                write(*,*) 'GenTurb: ', a,meshFaceCenters(1:3,iFace),ttime+dt
            end if
        end do
    end if
    
    do iFace = 1,nFaces
        Tsnew(iFace) = TemperaturePRho(Psnew(iFace),Rhosnew(iFace))
    end do
end subroutine CalcFacesVariables
